Waterwater Communities Manuscript Analysis

Supplemental Figure 1 - Coverage Analysis

SC- Hill number based coverage estimate, extrapolated rarefaction curve, considering the singleton to tentons (not sure this is a real word)? Coverage is the estimated OTU coverage calculated by dividing the number of observed OTUs by the Chao1 diversity estimate.
Chao1, is based on singletons and likely overestimates species diversity for 16S data. Therefore Chao1 based coverage estimates are conservative. Filtering samples with less than 100 sequences.

div_est <- plantObj  %>% MRcounts()  %>% DataInfo()
div_est$StageID <- factor(count_df$StageID)
p <- ggplot(div_est, aes(x = n, y = SC)) +
      geom_point(aes(color = StageID)) + scale_x_log10() + 
      geom_vline(aes(xintercept = 100), linetype = 2, color = "grey 60") + 
      theme_bw() +
      labs(x = "Number of Sequences",y = "Coverage") +
      theme(legend.position = "bottom") 
ggMarginal(p, type = "histogram", margins = "x")

Number of observed sequences compared to the estimated coverage. The vertical dotted line indicates the cut off for sample excluded from the analysis. Histogram at the top show the distribution of samples relative to the number of sequences per sample.

Excluding samples with less than 100 counts.

plantObj <- which(colSums(plantObj) > 100) %>% {plantObj[, .]} 
## creating phyloseq object
plant_phy <-plantObj %>% MRexperiment2biom(norm = TRUE,log = TRUE) %>% 
    import_biom2()
colnames(plant_phy@tax_table) <- c("taxa_OTU","Taxonomy","Kingdom",
                                   "Phylum","Class","Order",
                                   "Family","Genus","Species")
plant_phy@tax_table <- plant_phy@tax_table[,-2]

Supplemental Figure 2 - Alpha diversity

ww_div_df %>% filter(Stage == "Influent") %>% 
    ggplot() + geom_point(aes(x = plant_name, y = metric, color = plant_name)) + 
        facet_wrap(~diversity, scale = "free_y") + 
        theme_bw() + theme(axis.text.x = element_blank(), legend.position = "bottom") +
        labs(x = "Wastewater Treatment Plant", y = "Diversity Metric", color = "")

ww_inf_df <- ww_div_df %>% filter(Stage == "Influent")

fit_aov <- function(df) aov(metric~WW_TP, data = df)
test_kruskal <- function(df) kruskal.test(df$metric,factor(df$WW_TP))


inf_div_test <- ww_inf_df %>% group_by(diversity) %>% nest() %>% 
    mutate(aov_mod = map(data, fit_aov), krusk = map(data, test_kruskal))

aov_summary <- inf_div_test$aov_mod %>% map_df(glance) %>% 
    rename(f_stat = statistic, f_p = p.value) 
aov_summary$diversity <- inf_div_test$diversity

krusk_summary <- inf_div_test$krusk %>% map_df(glance) %>% 
    rename(kruskal = statistic, kruk_p = p.value) %>% 
    select(-parameter, -method)
krusk_summary$diversity <- inf_div_test$diversity

inf_div_test_tbl <- aov_summary %>% 
    select(diversity, logLik, AIC,BIC, f_stat, f_p) %>% left_join(krusk_summary)

inf_div_test_tbl %>% kable()

div_metrics <- c("shannon","simpson","specnumber")
for(i in 1:3){
    print(autoplot(inf_div_test$aov_mod[[i]]) + 
              ggtitle(paste0(div_metrics[i])))
}

Figure 1 - Relative Abundance

mrexp_inf <- plantObj[, which(pData(plantObj)$Stage_name %in% 
                                  c("Raw influent", 
                                    "Influent post screening"))]
mrexp_inf <- aggregateByTaxonomy(mrexp_inf, 
                                 lvl = "genus", 
                                 alternate = TRUE, 
                                 norm = TRUE, log = TRUE)

pData(mrexp_inf)$plant_name <- pData(mrexp_inf)$plant_name %>% factor()

mrexp_inf <- cumNorm(mrexp_inf, p = 0.75)
mrexp_inf <- filterData(mrexp_inf, present = 7, depth = 1)
# mrexp_inf <- cumNorm(mrexp_inf, p = 0.75)
s <- normFactors(mrexp_inf)
pd <- pData(mrexp_inf)

settings <- zigControl(maxit = 1, verbose = FALSE)
mod <- model.matrix(~plant_name, data = pd)
colnames(mod) <- levels(pd$plant_name)
res = fitZig(obj <- mrexp_inf, mod = mod, control = settings)
zigFit = res$fit
fit2 = eBayes(zigFit)
top_tbl <- topTableF(fit2,p.value = 0.075) 
top_feat <- top_tbl %>% rownames()

count_tbl <- MRcounts(mrexp_inf)
count_tbl <- count_tbl[rownames(count_tbl) %in% top_feat,]
count_df <- count_tbl %>% as.data.frame() %>% 
    rownames_to_column(var = "genus") %>% 
    gather("Seq_ID", "Count",-genus) %>% left_join(pd)

Differentially abundance genus between treatment plant influent samples. Might want to exclude MW_TP2 since there is only one sample.

kable(top_tbl)

count_df %>% 
    # mutate(plant_name = str_replace(plant_name, "Mid-Atlantic ", "MA-"),
    #        plant_name = str_replace(plant_name, "Midwest ", "MW-")) %>% 
ggplot() + 
    geom_point(aes(x = plant_name, y = Count, color = plant_name)) + 
    geom_line(aes(x = plant_name, y = Count, color = plant_name)) + 
    facet_wrap(~genus,nrow = 2) +
    theme_bw() + scale_y_log10() + labs(y = "Relative Abundance", 
                                        color = "Treatment Plant",
                                        x = "Treatment Plant") + 
    theme(legend.position = "bottom", 
          axis.text.x = element_text(angle = 90),
          strip.text =element_text(face = "italic"))

Supplemental Figure 3 - Alpha diversity

ww_div_inf_eff <- ww_div_df %>% filter(Stage %in% c("Influent","Effluent")) %>% 
    group_by(plant_name, Stage, Date, diversity) %>% 
    summarise(metric = median(metric))

ww_div_inf_eff %>% ungroup() %>% 
    mutate(Stage = factor(Stage, levels = c("Influent","Effluent")),
           Date = paste(plant_name, Date)) %>% 
    ggplot(aes(x = Stage, y = metric)) +
        geom_path(aes(group = Date), color = "grey60") +
        geom_point(aes(x = Stage, y = metric, color = plant_name)) +
        facet_wrap(~diversity, scale = "free_y") +
        labs(y = "Diversity Metric", color = "WWTP") + theme_bw() +
      theme(legend.position = "bottom")

Testing for differences

Using a paired t-test to compared the alpha diversity between influent and effluent samples, specnumber (number of observed OTUs) and shannon’s diversity index were significantly higher for influent compared to effluent. Shannon’s diversity index is close to the significance threshold of 0.05, and when accounting for multiple comparisons Shannon’s diversity index is no longer significant. The number of observed OTUs is biased by sequencing depth, as the effluent samples consistently have fewer reads than influent samples the observed difference is likely at least partially due to differences in sequencing depth.

inf_eff_paired <- ww_div_inf_eff %>% group_by(plant_name, Date, Stage, diversity) %>% 
    summarize(metric = median(metric)) %>% spread(Stage, metric) %>%
    filter(!is.na(Influent), !is.na(Effluent))

inf_v_ef_test_df <- data_frame()
for(i in c("shannon","simpson","specnumber")){
    inf_v_ef_test_df <- inf_eff_paired %>% 
        filter(diversity == i) %>% 
        {t.test(.$Influent, .$Effluent, paired = TRUE, alternative = "greater")} %>% 
        tidy() %>% mutate(diversity = i) %>% bind_rows(inf_v_ef_test_df)
}
inf_v_ef_test_df %>% select(diversity, estimate, statistic, p.value, parameter, conf.low, conf.high) %>% kable()

Relative Abundance

mrexp_inf_eff <- plantObj[, which(pData(plantObj)$Stage_name %in% 
                                  c("Raw influent", 
                                    "Influent post screening","Effluent"))]
mrexp_inf_eff <- aggregateByTaxonomy(mrexp_inf_eff, 
                                 lvl = "genus", 
                                 alternate = TRUE, 
                                 norm = TRUE, log = TRUE)

pData(mrexp_inf_eff)$plant_name <- pData(mrexp_inf_eff)$plant_name %>% factor()
pData(mrexp_inf_eff)$Stage_name <- pData(mrexp_inf_eff)$Stage_name %>% 
      ifelse(. %in% c("Raw influent","Influent post screening"),"Influent",.) %>% 
      factor(levels = c("Influent","Effluent"))

mrexp_inf_eff <- cumNorm(mrexp_inf_eff, p = 0.75)
mrexp_inf_eff <- filterData(mrexp_inf_eff, present = 7, depth = 1)
# mrexp_inf_eff <- cumNorm(mrexp_inf_eff, p = 0.75)
s <- normFactors(mrexp_inf_eff)
pd <- pData(mrexp_inf_eff)

settings = zigControl(maxit = 1, verbose = FALSE)
mod = model.matrix(~Stage_name, data = pd)
colnames(mod) = levels(pd$Stage_name)
res = fitZig(obj = mrexp_inf_eff, mod = mod, control = settings)
zigFit = res$fit
fit2 <- eBayes(zigFit)
top_tbl <- topTableF(fit2,p.value = 0.075) 
top_feat <- top_tbl %>% rownames()

count_tbl <- MRcounts(mrexp_inf_eff)
count_tbl <- count_tbl[rownames(count_tbl) %in% top_feat,]
count_df <- count_tbl %>% as.data.frame() %>% 
    rownames_to_column(var = "genus") %>% 
    gather("Seq_ID", "Count",-genus) %>% left_join(pd) %>% 
    mutate(plant_date = paste0(Date, plant_name))

kable(top_tbl)

Figure 2

ggplot(count_df) + 
    geom_path(aes(x = Stage_name, y = Count, group = plant_date), 
              color = "grey60") + 
    geom_point(aes(x = Stage_name, y = Count, color = plant_name)) + 
    facet_wrap(~genus, nrow = 2) +
    theme_bw() + scale_y_log10() + 
    labs(y = "Relative Abundance", x = "Stage Name", color = "Treatment Plant") +
    theme(legend.position = "bottom", 
          strip.text = element_text(face = "italic"))

Treatment Process Comparison

pl <- subset_samples(plant_phy, study_type == "plant")

pl_pcoa_list <- list()
pl_nmds_list <- list()
for (i in dist_methods) {
    # Calculate distance matrix
    iDist <- distance(pl, method = i)
    # Calculate ordination
    pl_pcoa_list[[i]] <- ordinate(pl, method = "PCoA", distance = iDist)
    pl_nmds_list[[i]] <- ordinate(pl, method = "NMDS", distance = iDist)
}

Figure 3 - Ordination Plots

# plot_ordination(pl, pl_pcoa_list[[1]], color= "Stage", title = "Bray Curtis PCoA") + theme_bw() + stat_ellipse()
# plot_ordination(pl, pl_nmds_list[[1]], color= "Stage", title = "Bray Curtis NMDS") + theme_bw() + stat_ellipse()
# stressplot(pl_nmds_list[[1]])

# plot_ordination(pl, pl_pcoa_list[[2]], color= "Stage", title = "Jaccard PCoA") + theme_bw()
# plot_ordination(pl, pl_nmds_list[[2]], color= "Stage", title = "Jaccard NMDS") + theme_bw()
# stressplot(pl_nmds_list[[2]])

Testing for Differences

iDist <- distance(pl, method = "bray")
anosim(iDist,grouping = factor(pl@sam_data$Stage))

## 
## Call:
## anosim(dat = iDist, grouping = factor(pl@sam_data$Stage)) 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.5632 
##       Significance: 0.001 
## 
## Permutation: free
## Number of permutations: 999

Treatment plant stages are significantly different testing for pairwise differences

## betadisper has a post-hoc test
iDist <- distance(pl, method = "bray")
beta_fit <- betadisper(iDist,group = factor(pl@sam_data$Stage))

anova(beta_fit) %>% tidy() %>% kable()

TukeyHSD(beta_fit)$group %>% as.data.frame() %>% 
    rownames_to_column(var = "comparison") %>% 
    filter(`p adj` < 0.075) %>% kable(digits = 3)

Spray Irrigation Stage Comparison

si <- subset_samples(plant_phy, study_type == "spray")
si_pcoa_list <- list()
si_nmds_list <- list()
for (i in dist_methods) {
    print(i)
    # Calculate distance matrix
    iDist <- distance(si, method = i)
    # Calculate ordination
    si_pcoa_list[[i]] <- ordinate(si, method = "PCoA", distance = iDist)
    si_nmds_list[[i]] <- ordinate(si, method = "NMDS", distance = iDist)
}

## [1] "bray"
## Run 0 stress 0.08628405 
## Run 1 stress 0.08628405 
## ... Procrustes: rmse 9.610914e-07  max resid 2.007411e-06 
## ... Similar to previous best
## Run 2 stress 0.3829659 
## Run 3 stress 0.2558146 
## Run 4 stress 0.08628405 
## ... Procrustes: rmse 2.892117e-06  max resid 9.106018e-06 
## ... Similar to previous best
## Run 5 stress 0.08628405 
## ... New best solution
## ... Procrustes: rmse 1.062147e-06  max resid 1.871165e-06 
## ... Similar to previous best
## Run 6 stress 0.09761258 
## Run 7 stress 0.1074574 
## Run 8 stress 0.1601557 
## Run 9 stress 0.09761272 
## Run 10 stress 0.08628405 
## ... Procrustes: rmse 2.545378e-06  max resid 4.930609e-06 
## ... Similar to previous best
## Run 11 stress 0.0964332 
## Run 12 stress 0.1600244 
## Run 13 stress 0.1074574 
## Run 14 stress 0.1663974 
## Run 15 stress 0.09643206 
## Run 16 stress 0.08628405 
## ... Procrustes: rmse 8.360176e-07  max resid 1.547145e-06 
## ... Similar to previous best
## Run 17 stress 0.1609997 
## Run 18 stress 0.1074574 
## Run 19 stress 0.09636965 
## Run 20 stress 0.0976128 
## *** Solution reached
## [1] "jaccard"
## Run 0 stress 0.08628405 
## Run 1 stress 0.08628405 
## ... New best solution
## ... Procrustes: rmse 1.580754e-06  max resid 5.108713e-06 
## ... Similar to previous best
## Run 2 stress 0.09761348 
## Run 3 stress 0.3856545 
## Run 4 stress 0.08628405 
## ... Procrustes: rmse 3.412486e-06  max resid 1.333437e-05 
## ... Similar to previous best
## Run 5 stress 0.09761607 
## Run 6 stress 0.1074674 
## Run 7 stress 0.08628405 
## ... Procrustes: rmse 2.537932e-06  max resid 6.172675e-06 
## ... Similar to previous best
## Run 8 stress 0.1601557 
## Run 9 stress 0.09761644 
## Run 10 stress 0.09628348 
## Run 11 stress 0.1617694 
## Run 12 stress 0.08628405 
## ... Procrustes: rmse 2.067053e-06  max resid 5.999207e-06 
## ... Similar to previous best
## Run 13 stress 0.08628405 
## ... Procrustes: rmse 1.855991e-06  max resid 6.65927e-06 
## ... Similar to previous best
## Run 14 stress 0.09636983 
## Run 15 stress 0.08628405 
## ... Procrustes: rmse 4.220688e-06  max resid 1.290558e-05 
## ... Similar to previous best
## Run 16 stress 0.1601557 
## Run 17 stress 0.08628405 
## ... New best solution
## ... Procrustes: rmse 6.392372e-07  max resid 2.221082e-06 
## ... Similar to previous best
## Run 18 stress 0.08628405 
## ... Procrustes: rmse 1.161473e-06  max resid 3.827379e-06 
## ... Similar to previous best
## Run 19 stress 0.09625276 
## Run 20 stress 0.08628405 
## ... Procrustes: rmse 1.961568e-06  max resid 4.142391e-06 
## ... Similar to previous best
## *** Solution reached

Ordination Plots

#plot_ordination(si, si_pcoa_list[[1]], color= "Stage", title = "Bray Curtis PCoA") + theme_bw()
#plot_ordination(si, si_nmds_list[[1]], color= "Stage", title = "Bray Curtis NMDS") + theme_bw()
stressplot(si_nmds_list[[1]])

#plot_ordination(si, si_pcoa_list[[2]], color= "Stage", title = "Jaccard PCoA") + theme_bw()
#plot_ordination(si, si_nmds_list[[2]], color= "Stage", title = "Jaccard NMDS") + theme_bw()
stressplot(si_nmds_list[[2]])

Testing for Differences by Stage

Spray irrigation stages are significantly different.

## betadisper has a post-hoc test
iDist <- distance(pl, method = "bray")
beta_fit <- betadisper(iDist,group = factor(pl@sam_data$Stage))

anova(beta_fit) %>% tidy() %>% kable()

TukeyHSD(beta_fit)$group %>% as.data.frame() %>% 
    rownames_to_column(var = "comparison") %>% 
    filter(`p adj` < 0.075) %>% kable(digits = 3)

Supplemental Figure 4 - Alpha diversity

ww_div_treat <- ww_div_df %>% filter(study_type == "plant") %>% 
    group_by(plant_name, Stage, StageID, Date, diversity) %>% 
    summarise(metric = median(metric))

ww_div_treat %>% ungroup() %>% 
    mutate(Stage = factor(Stage,levels = c("Influent","ActivatedSludgeReactor", 
                                           "PostAeration","SecondaryClarifier",
                                           "CellB","Effluent")), 
           Date = paste(plant_name, Date)) %>% 
    ggplot(aes(x = Stage, y = metric)) +
        geom_line(aes(group = Date), color = "grey60") +
        geom_point(aes(x = Stage, y = metric)) +
        facet_grid(diversity~plant_name, scale = "free") +
        labs(y = "Diversity Metric", color = "WWTP") + theme_bw() +
        theme(axis.text.x = element_text(angle = 90))

fit_aov <- function(df) aov(metric~Stage, data = df)
treat_div_test <- ww_div_treat %>% filter(plant_name == "Mid-Atlantic WWTP1") %>% 
    group_by(diversity) %>% nest() %>% mutate(aov_mod = map(data, fit_aov))

aov_summary <- treat_div_test$aov_mod %>% map_df(glance) %>% 
    rename(f_stat = statistic, f_p = p.value) 
aov_summary$diversity <- inf_div_test$diversity

aov_summary %>% select(diversity, logLik, AIC,BIC, f_stat, f_p) %>% kable()

div_metrics <- c("shannon","simpson","specnumber")
for(i in 1:3){
    print(autoplot(treat_div_test$aov_mod[[i]]) + 
              ggtitle(paste0(div_metrics[i])))
}

Relative Abundance

mrexp <- aggregateByTaxonomy(plantObj, 
                                 lvl = "genus", 
                                 alternate = TRUE, 
                                 norm = TRUE, log = TRUE)
pData(mrexp)$plant_name <- pData(mrexp)$plant_name %>% factor()
pData(mrexp)$Stage_name <- pData(mrexp)$Stage_name %>% 
      ifelse(. %in% c("Raw influent","Influent post screening"),"Influent",.) %>% 
      factor(levels =c("Influent", "Activated Sludge Reactor", "Post aeration", 
                       "Secondary clarifier", "Cell B","Effluent"))

plant_zigFit <- function(mrexp, plant){
    mrexp_plant <- mrexp[,which(pData(mrexp)$plant_name == plant)]
    pData(mrexp_plant)$Stage_name <- factor(pData(mrexp_plant)$Stage_name)
    mrexp_plant <- cumNorm(mrexp_plant, p = 0.75)
    mrexp_plant <- filterData(mrexp_plant, 
                              present = floor(dims(mrexp_plant)[2]/2), 
                              depth = 1)
    # mrexp_plant <- cumNorm(mrexp_plant, p = 0.75)
    pd <- pData(mrexp_plant)
    
    settings = zigControl(maxit = 1, verbose = FALSE)
    mod = model.matrix(~Stage_name, data = pd)
    colnames(mod) = levels(pd$Stage_name)

    res = fitZig(obj = mrexp_plant, mod = mod, control = settings)
    zigFit = res$fit
    finalMod = res$fit$design
    fit2 = eBayes(zigFit)
    top_tbl <- topTableF(fit2,p.value = 0.075) 
    top_feat <- top_tbl %>% rownames()

    count_tbl <- MRcounts(mrexp_plant)
    count_tbl <- count_tbl[rownames(count_tbl) %in% top_feat,]
    count_df <- count_tbl %>% as.data.frame() %>%
        rownames_to_column(var = "genus") %>%
        gather("Seq_ID", "Count",-genus) %>% left_join(pd)

   list(count_df = count_df, top_tbl = top_tbl)
}

Figure 4-7

make_ww_tp_plot <- function(plant_zig){
    ggplot(plant_zig$count_df) +
                geom_point(aes(x = Stage_name, y = Count, color = Stage_name)) +
                facet_wrap(~genus, nrow = 1) +
                theme_bw() + scale_y_log10() +
                labs(y = "Rel. Abu.", color = "Treatment Process") +
                theme_bw() +
              theme(axis.text.x = element_blank(),
                    axis.title.x = element_blank(),
                    strip.text = element_text(face = "italic"))
}

ww_tp <- pData(mrexp)$plant_name %>% 
    .[!grepl("SI", .)] %>% unique() %>% as.character()

ww_tp

## [1] "Mid-Atlantic WWTP1" "Midwest WWTP1"      "Midwest WWTP2"     
## [4] "Mid-Atlantic WWTP2"

ww_tp_plots <- ww_tp %>% set_names(ww_tp) %>% 
    map(plant_zigFit, mrexp = mrexp) %>% map(make_ww_tp_plot)

ggarrange(ww_tp_plots[[1]], ww_tp_plots[[4]], 
          ww_tp_plots[[2]], ww_tp_plots[[3]],
          labels = c("A", "B", "C", "D"),
          ncol = 1, nrow = 4,legend = "bottom")

ggsave("diff_abu_treat_process.tiff", width = 10, height = 8, dpi = 300)

for(i in ww_tp){
    plant_zig <- plant_zigFit(mrexp, i)
    print(paste("##### ",i))
    print(kable(plant_zig$top_tbl))
}

[1] “##### Mid-Atlantic WWTP1”

Supplemental Figure S5 - Alpha diversity

spray_div <- ww_div_df %>% filter( WW_TP %in% c("MA_TP1","MA_SI1")) %>% 
    filter(Stage %in% c("Influent","Effluent", "BeforeUVTreatment", 
                        "AfterUVTreatment","HoldingPondInlet", 
                        "PumpHouseInlet")) %>% 
    group_by(plant_name, Stage, StageID, Date, diversity) %>% 
    summarise(metric = median(metric))

spray_div %>% ungroup() %>% 
    mutate(Stage = factor(Stage,levels = c("Influent","Effluent","BeforeUVTreatment", 
                        "AfterUVTreatment","HoldingPondInlet", 
                        "PumpHouseInlet"))) %>% 
    ggplot(aes(x = Stage, y = metric)) +
        geom_line(aes(group = Date), color = "grey60") +
        geom_point(aes(x = Stage, y = metric)) +
        facet_grid(diversity~., scale = "free") +
        labs(y = "Diversity Metric", color = "WWTP") + theme_bw() +
        theme(axis.text.x = element_text(angle = 90))

fit_aov <- function(df) aov(metric~Stage, data = df)
spray_div_test <- spray_div %>% group_by(diversity) %>% 
    nest() %>% mutate(aov_mod = map(data, fit_aov))

aov_summary <- spray_div_test$aov_mod %>% map_df(glance) %>% 
    rename(f_stat = statistic, f_p = p.value) 
aov_summary$diversity <- spray_div_test$diversity

aov_summary %>% select(diversity, logLik, AIC,BIC, f_stat, f_p) %>% kable()

div_metrics <- c("shannon","simpson","specnumber")
for(i in 1:3){
    print(autoplot(spray_div_test$aov_mod[[i]]) + 
              ggtitle(paste0(div_metrics[i])))
}

for(i in c("shannon", "specnumber")){
    print(i)
    tukey <- spray_div %>% filter(diversity == i) %>% {aov(metric~Stage, data = .)} %>% TukeyHSD()
    tidy(tukey) %>% filter(adj.p.value < 0.1) %>% print()
}

## [1] "shannon"
##    term                comparison estimate   conf.low conf.high
## 1 Stage Influent-HoldingPondInlet 1.374929 0.06097935   2.68888
##   adj.p.value
## 1  0.03698308
## [1] "specnumber"
##    term                 comparison estimate conf.low conf.high adj.p.value
## 1 Stage  Influent-AfterUVTreatment 1271.125 287.8470  2254.403 0.006898620
## 2 Stage Influent-BeforeUVTreatment 1161.500 178.2220  2144.778 0.014936608
## 3 Stage  Influent-HoldingPondInlet 1384.083 486.4774  2281.689 0.001193148

Relative Abundance

mrexp_spray <- plantObj[,which(pData(plantObj)$WW_TP %in% c("MA_TP1","MA_SI1") & 
                              pData(plantObj)$Stage %in% c("RawInfluent","PostScreeningInfluent","Effluent",
                                                            "BeforeUVTreatment", "AfterUVTreatment",
                                                            "HoldingPondInlet","PumpHouseInlet"))]

pData(mrexp_spray)$Stage_name <- pData(mrexp_spray)$Stage %>% 
      ifelse(. %in% c("RawInfluent","PostScreeningInfluent"),"Influent",.) %>% 
      factor(levels =c("Influent","Effluent", "BeforeUVTreatment", "AfterUVTreatment",
                       "HoldingPondInlet","PumpHouseInlet"))

mrexp_spray <- aggregateByTaxonomy(mrexp_spray, 
                                 lvl = "genus", 
                                 alternate = TRUE, 
                                 norm = TRUE, log = TRUE)

mrexp_spray <- cumNorm(mrexp_spray, p = 0.75)
mrexp_spray <- filterData(mrexp_spray, present = 7, depth = 1)
# mrexp_spray <- cumNorm(mrexp_spray, p = 0.75)
s <- normFactors(mrexp_spray)
pd <- pData(mrexp_spray)

settings = zigControl(maxit = 1, verbose = FALSE)
mod = model.matrix(~Stage_name, data = pd)
colnames(mod) = levels(pd$Stage_name)
res = fitZig(obj = mrexp_spray, mod = mod, control = settings)
zigFit = res$fit
fit2 <- eBayes(zigFit)
top_tbl <- topTableF(fit2,p.value = 0.075) 
top_feat <- top_tbl %>% rownames()

count_tbl <- MRcounts(mrexp_spray)
count_tbl <- count_tbl[rownames(count_tbl) %in% top_feat,]
count_df <- count_tbl %>% as.data.frame() %>% 
    rownames_to_column(var = "genus") %>% 
    gather("Seq_ID", "Count",-genus) %>% left_join(pd) %>% 
    mutate(Stage = factor(Stage,levels = c("RawInfluent","Effluent","BeforeUVTreatment", 
                    "AfterUVTreatment","HoldingPondInlet", 
                    "PumpHouseInlet")))

kable(top_tbl)

ggplot(count_df) + 
    geom_line(aes(x = Stage, y = Count, group = Date), color = "grey60") + 
    geom_point(aes(x = Stage, y = Count)) + 
    facet_wrap(~genus, nrow = 2) +
    theme_bw() + scale_y_log10() + labs(y = "Relative Abundance") +
    theme(axis.text.x = element_text(angle = 90))